Numerical study of quasiparticle lifetime in quantum dots 
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The decay rate of quasiparticles in quantum dots is studied through the real time calculation of 
the single-particle Green function in the self-consistent approximation. The method avoids exact 
diagonalization, transforming the problem into a system of coupled non-linear integral equations 
which may be solved iteratively. That allows us to study systems larger than previously treated 
in the literature. Our results for the inverse participation ratio of many-body states show that the 
threshold energy for the quasiparticle disintegration is E* ~ \f§A. The derealization transition is 
soft rather than sharp. Three different regimes as function of the effective interaction strength may 
be clearly identified at high energies. 



I. INTRODUCTION 

Isolated quantum dots provide a unique opportunity 
to study the interplay between phase coherence, disor- 
der, and electron-electron interaction effects in a confined 
geometry. These systems are realized experimentally 
in several forms; the best known example are the elec- 
tron islands formed by the electrostatic lateral confine- 
ment of a two-diaiensional electron gas in GaAs / GaAlAs 
heterostructuresJd Among the important issues recently 
investigated is the lifetime of quasiparticle excitations 
that can be associated with the broadening of tunnel- 
ing conductance peaks. In experiment reported in Ref. 
|2|, for instance, few narrow peaks were observed above 
the Fermi energy, followed by a broad continuum. This 
result has initiated a theoretical discussion about the na- 
ture of the excited states in confined (but large) fermionic 
systems. 

In this paper we study numerically the problem of 
quasiparticle decay rate in finite fermionic systems with 
chaotic single-particle dynamics. In most quantum dots, 
the presence of either disorder or complex boundaries 
lead to single-particle eigenstates displaying universal 
statistical properties. The universality takes place for 
states within the energy range of order gA near the 
Fermi energy, where g 1 is the dimensionless con- 
ductance of the system (i.e, conductance in units e 2 /K) 
and A is the mean spacing between neighboring single- 
particle eigenenergies. The characteristic features of such 
states are level repulsion and spectral rigidityJa Moreover, 
the two-body (short-ranged) interaction matrix elements 
taken in the single-particle eigenbasis fluctuate according 
to Gaussian distribution whose standard deviation scales 
as 1/ga 

Ref. [5] proposed to map the problem of quasiparticle 
lifetime onto an Anderson problem in the many-body 
Fock space. Upon such mapping, many-body states of 



the noninteracting system are considered as "sites" of 
the Anderson lattice, while the interactions provide hop- 
ping between them. Since both the "on-site" energies 
and the hopping matrix elements are random (disorder 
dependent), one may expect a localization-delocalization 
transition to take place for certain values of the system 
parameters. For excitations with energy e lower than 
the threshold energy E* , the exact many-body states 
are "localized" in the single-particle basis. As a result, 
an excited quasiparticle evolves into a superposition of 
only few low-energy many-body eigenstates. The spec- 
tral function consists of a finite number of 5-functional 
peaks and the quasiparticle lifetime is essentially infinite 
in this regime. At the energies above the derealization 
threshold, e > E* , a large number of basis single-particle 
states contribute to a given exact many-body states. This 
eventually leads to spectral functions of the Breit-Wigncr 
form whose width (inverse lifetime) may be evaluated em- 
ploying the golden rule (GR). 

To work out details of this crossover, Ref. || employed 
an approximate mapping of the many-body space of a 
fermionic system onto a Cayley tree (CT). Each genera- 
tion of the CT corresponds to a set of many-body states 
with a given number of electron- hole pairs. Interactions 
(hopping) provide coupling between generations by excit- 
ing new electron-hole pairs. The advantage of such anal- 
ogy is that the Anderson localization problem on the CT 
may be solved exactly.l§El Such solution indicates t he exis- 
tence of a sharp crossover at the energy E* ~ Ay/gjhig , 
where the participation ratio P jumps abruptly from a 
value close to unity in the localized phase (e < E*) to 
infinity in the delocalized phase (e > E*). The other im- 
portant feature of the exact CT solution is the existence 
of a second characteristic energy E** ~ A^/g > E* . For 
excitation energies e > E**, the spectral density is of the 
Breit-Wigner form and the GR applies. On the other 
hand, for energies E* < e < E** , the spectral density 
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consists of a very large number of narrow, but separate 
peaks (though the participation ratio is infinite) . 

The disadvantage of the CT mapping is that it is by no 
means exact. Indeed, it disregards potentially important 
effects such as existence of loops in the many-body space, 
correlations between on-site energies of different genera- 
tions, and the fact that the actual number of accessible 
many-body states at a given energy is finite. Therefore, 
the predictions of the exact CT solution should be taken 
with a grain of salt. In fact, it is clear that no abrupt 
transition may happen, since the number of many-body 
states at the energy close to the critical one is exponen- 
tially large, but finite. As a result, the very existence 
of the transition, its quantitative features, and the char- 
acteristic energy scales are not firmly established bypth^i 
CT picture. Several theoretical and numerical worksQ~t3 
were devotad to clarification of these issues. One numer- 
ical studyO claimed that for small systems the transi- 
tion is absent and the quasiparticlc lifetime decay can be 
well described by the GR. That conclusion does not ap- 
ply to larger systems, where a breakdown of the GR is 
observedlij. In large systems, one can explore a broader 
range of values of g, as well as higher excitation energies. 
Nevertheless, no signal of an abrupt transition was found. 
Essentially, all numerical work so far has relied on the 
exact diagonalization of the full many-body Hamiltonian 
(thus limited size Fock space), although under different 
assumptions regarding the structure of the two-body ma- 
trix elements. 

Here we employ an alternative approach suggested by 
Levitovii3 The method avoids matrix inversion or di- 
agonalization and allows the investigation of large size 
systems. It is based on the calculation of the single- 
particle Green's function in real time. Adopting the 
self-consistent approximation for the self-energy, one can 
transform the problem into a non-linear integral equa- 
tion which may be solved iteratively. The inverse partic- 
ipation ratio that measures the overlap between a one- 
particle state and the many-body eigenstates may be ob- 
tained directly from the long-time behavior of the Green's 
function. Our simulations involved up to 50 particles and 
150 basis states (thus with filling factor 1/3; we deal with 
the spinless fermions for simplicity). The results confirm 
conclusions of Ref. [13], namely, that there is no sharp 
transition of the participation ratio (contrary to the case 
of an ideal CT). However, our estimate for the thresh- 
old energy at which quasiparticle states delocalize does 
not match the prediction of Ref. 0. Instead, our results 
are consistent with the scaling of the threshold energy 
as A^/g. We see no indications of other scalings such as 
those predicted in Ref. |lj. 

At the same time, our numerical data is consistent with 
the existence of the three different energy regimes, in 
agreement with Ref. [| At small energies, single-particle 
states overlap with only a finite number of many-body 
states and quasiparticles do not decay. In the the oppo- 
site limit, of large energy, the quasiparticles decay rapidly 
into the continuum of many-body states. The decay rate 



follows closely the GR prediction. In the intermediate 
energy range, the quasiparticle is neither well defined 
(Fock-space localized), nor short lived, and GR does not 
apply. 

This paper is organized as follows. In Sec. || we for- 
mulate the method used in the numerical simulations. 
The results are presented in Sec. [II. Discussion and 



conclusions are left to Sec. IV 



II. REAL-TIME GREEN'S FUNCTION 
APPROACH 

The Hamiltonian of an isolated quantum dot with N 
spinless electrons may be separated into two parts, 



H — Ho + Hi , 



(1) 



where Hq contains one-body terms, such as the kinetic 
energy and the background potential, and Hi represents 
the two-body interactions. We denote £j and <f>j(x) the 
single-particle eigenvalues and eigenfunctions of H , re- 
spectively. In this eigenbasis, one may write 



£j CjCj 



and 



Hi = 



j<k,m<l 



v jk c m c l c J c k> 



(2) 



(3) 



where the fermionic operators Cj and cj are, respectively, 
annihilation and creation operators of an electron on the 
state j. The interaction matrix elements are given by 



VjF = J dx J dx' <f>* m {x)4>i{x') V(x - x') <j> k {x)^{x'). 

(4) 

Below, we consider the case of short-range interactions 
of the form 



V(x) — XAL d S^(x) 



(5) 



where L d is the system volume and A is the dimension- 
less interaction strength. Diagonal and off-diagonal ma- 
trix elements of the interaction have distinct statistical 
propertiesJj The diagonal ones, 



V] k k = V%=\AL« ldx\MxW\MxW, (6) 



have non-zero average value, Vj k ~ AA, while fluctua- 
tions around this mean value are small as g . On the 
other hand the off-diagonal matrix elements 

Vjr — AA L d j dx cf>* m (x)<fi {x)<t> k (x)^ (x) (7) 
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are zero in average, VjjJ 1 = 0, while their typical value is 
\Vj™\ = XA/g for any choice of different j, k, I, to within 
the universal range Ag.S 

In the absence of Hi , the ground state of the system is 
formed by the Slater determinant of the N lowest energy 
single particle eigenstates. That is, for eigenvalues ar- 
ranged in ascending order, all states below and including 
the Fermi energy £p = sn are occupied and the remain- 
ing are empty. Excitations are formed by promoting one 
or more particles from the Fermi sea to the empty states 
above Ef, leaving the holes behind. For instance, for an 
excited cigcnstatc formed by to particle-hole pairs, 



b'2m, ■ • ■ ,jm+l]jm> ■ ■ ■ j Jl)o 



ct c f 



xc jm --- Cjl \FS) . (8) 



(Here, the subscript is used to stress the absence of 
interactions.) 



A. The single—particle Green's function 

Let us consider the (zero-temperature) time-ordered, 
single-particle Green's function (Ti = 1): 

G^t) = -i(FS|r{ Ci (t)ct(0)}|FS), (9) 

where |FS) is the exact many-body ground state of the 
full Hamiltonian and T denotes the time ordering. 

It is useful to review briefly the non-interacting case. 
The Green's function can then be written in the following 
form, 



Gf\t) 



-i6(t) e-'*** , Ej > E F 

-ie-jt 



i6(-t)e 



ej < e F 



(10) 



where 9(t) is the Heaviside step function. It is straight- 
forward to show that 



(id t -Ej)Gf\t) = S(t), 



.(0) 



with the initial conditions 

gJV) = -», 

Gf } (0-) = i, 



After the Fourier transform of Eq. Oil), one finds 



e- <e (0) 

fcj ^ t F 



G?\E) 



1 



E-£i 



i0 + sgn(ej) ' 



(11) 

(12) 
(13) 

(14) 



Thus, Gj(E) has a pole whose real part is exactly the 
j-th single-particle eigenenergy. 

In the presence of interactions, Hi, the ground state is 
a superposition of a large number of Slater determinants. 
The same applies for the excited states, which are formed 



by a superposition of many particle-hole states. Never- 
theless, many aspects of the interactions, not related to 
the correlation effects, may be taken into account in the 
Hartree-Fock approximation. For instance, the diagonal 
part of Hi (related to matrix elements V? k and V*£) 
can be incorporated into Ho, creating new sets of single- 
particle eigenstates {4>j} and eigenenergies {e^ }, with 
the statistical properties similar to the noninteracting 
ones.ta For our purposes, it is unnecessary to distinguish 
this new set from the original one and the Hartree-Fock 
superscript is dropped. By the same token, one may 
drop the diagonal elements. The off-diagonal ones are 
kept and are solely responsible for the finite lifetime of 
the quasiparticles. 

The exact single-particle Green's function may be 
expanded perturbatively in ascending powers of A. 
Through standard diagrammatic techniques, one arrives 
at the Dyson equation which generalizes Eq. (|ll|) onto 
the interacting case, 

(id t - sj) Gj{t) = S(t) + J dt' E^(t - t') Gjit'). (15) 

The self-energy is the sum of all irreducible dia- 

grams and contains the information about the structure 
of the excited states. In the frequency representation, 
one obtains 



1 



E-Bi-^E) 



(16) 



Thus, the self-energy shifts the pole of the Green's func- 
tion with respect to the single-particle eigenenergy. 

In the scheme proposed in Ref. the self-energy in- 
volves only a limited class of diagrams. In order to moti- 
vate the approximation involved, let us view each excited 
many-body state formed by a Slater determinant of N 
eigenstates {<j>j } as a site in Fock space. The two-body in- 
teraction term acts as a hopping matrix elements between 
sites, while the energy at each site is formed by a sum of 
the energies of the single-particle eigenstates present in 
the corresponding Slater determinant. The sites are con- 
sidered as nearest neighbors when they differ only by one 
single-particle eigenstate. In fact, a more natural classi- 
fication of the excited manjj-body states uses the number 
of particle and hole states£3 Taking the noninteracting 
Fermi reference, we call a many-body state of 

class m (m > 0) when it is formed by m particle and 
to — 1 hole excitations. In Ref. |5|, hopping between sites 
of the same class were not allowed; a site in class m could 
connect to all accessible sites in class to + 1 , but to only 
one site in class to — 1 . It is not difficult to see that these 
restrictions lead to the topology of a CT. It is important 
to notice that the coordination number of the tree at a 
given site depends on the site class and the excitation 
energy; eventually, for a sufficiently high class number, 
one runs out of single-particle states and further branch- 
ing is interrupted. This fact was not properly taken into 
account in Ref. | whwi applying the results of Anderson 
localization problemJa As a result, while the Fock space 
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and the CT localization problems share a similar math- 
ematical structure, they are not fully equivalent. 

The same physics as in Ref. [|, with the proper branch- 
ing number and on-site energy correlations, is fully en- 
coded in the following self-consistent approximation for 
the self-energy, 



V kl 



kit. 



G k (t)Gi(t) G m (~t). 



(17) 



The resulting self-consistent equation for the Green's 
function is represented diagrammatically in Fig. (|l]). 
Equation (|l7j ) for the self-energy is the main approxi- 
mation adopted in the present paper. We believe that 
it encompasses the essential physics of the model. In 
particular, it includes the "correct" CT structure with 
floating coordination number, as well as some of the pos- 
sible loops in the many-body space. On the other hand, 
we have to stress that it is certainly not exact and some 
of the correlation effects are neglected by adopting Eq. 

©■ 



+ 



L = 




FIG. 1. Dyson equation and self-energy. 



B. The retardation ansatz 

For our purposes, it is usefulEl to work with the real- 
time formulation, as in Eqs. ( |l5| ) and ([l7]), instead of 
the energy representation. In order to allow for an it- 
erative solution of the non-linear coupled equations, we 
decompose the time-ordered Green's in particle-like and 
hole-like parts, depending on whether the index reference 
state j is above (sj > £f) or below (£j < Ef) the Fermi 
energy, respectively. That is, 



-ie Bit gj{t), sj > e F 
e-^fji-t), e 3 <s F 



(18) 



It is now straightforward to obtain a set of integro- 
differential equations for the functions gj(t) and fj(t). 
Let us, for convenience, assume ef = 0, such that j > 
corresponds to particle-like excitations, whereas j < 
corresponds to hole-like excitations. 

Using Eq. (Jisj) into Eqs. (||) and (0) for the case 
j > 0, we find that 



d 9j(t) 
dt 



= S(t)- / dt' S 3 {t-t') 9j (t') 



(19) 



with the kernel 



jrn 
kl 



Jt(ej -Ejj -ej+em) 



9k{t)gi{t)f m {t). 



k,l>0 
m<0 



Similarly, for j < we obtain 



dfjjt) 
dt 



5(t) - / dt' Rjlt-t^fjit') 



(20) 



(21) 



with the kernel 



RAt) = E \ v i 



jm 
kl 



^-it(sj-e k -ei+s„ 



k,l<0 
m>0 



] fk(i)fi{t)g m {t). 



(22) 



So far no additional approximation to Eq. ( [l7| ) have 
been made. However, the self-consistent equations are 
still too complicated to be treated numerically. To pro- 
ceed, we can explore the fact that the interactions are 
moderately weak (A ~ 1). Thus, the matrix elements 
Vfc/™ mix strongly only states that are nearly on shell, 
such that Se — e.j — Sk — £i + E m is of order of XA/g. 
Since we are interested in, say, particle-like states with 
the typical energy e ^> Se we can disregard their admix- 
ture with the hole-like states and vice-versa. Neglect- 
ing such mixing amounts to looking for the solutions of 
the coupled integro-differential equations in the retarded 
form. (As matter of fact we disregard the difference be- 
tween the retarded and time-ordered Green's functions 
for sufficiently large energies.) We thus set the initial 
conditions as stated above and substitute Eqs. ( |l9|) and 
©by 



and 



dg 3 (t) 

dt 



dm 

dt 



dt' Sj(t - t') gj {t') 



(23) 



[ dt'Rjit-t')^'), (24) 
Jo 



respectively, for t > 0. The computational advantage 
of the new equations over the original ones is obvious. 
For a given t, the integration runs over previous times 
only. Therefore, one can use the equations to iterate 
the initial condition forward in time to a certain given 
time t > 0. The error made in taking the retarded func- 
tions instead of the time-ordered one is measured by the 
parameter \V k r j n \/Ej. In the region of the expected tran- 
sition, Ej ~ Ay 7 !?, this parameter is small. We believe 
therefore that the retardation ansatz, Eqs. ( p3|) and 
(p4|), carries the same physical information as the self- 
consistent approximation, Eq. (p~7|) . 
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C. Inverse participation ratio and real-time 
quasiparticle decay 

Let us discuss characterization of the transition be- 
tween the extended and localized states in the Fock space 
by the real-time single-particle Green's function method. 
Going back to Eq. (w, for t > 0, one finds 



;X>M*)I«) (a c t(0) 



,it(E -E a ) 



(25) 



where \a) denotes an cigcnstatc of the full Hamiltonian, 
H, with corresponding eigenenergy E a . As a result: 

\GM 2 =^|(0| Ci |a)| 4 + 2^K0| Cj |a)| 2 |(0| Ci |^)| 2 

a a>/3 

x cos [{E a ~ E f3 )t] . (26) 

Upon ensemble averaging, the oscillating term disap- 
pears. Indeed, since E a — Ep behaves as a random vari- 
able, at long enough times the average of the oscillating 
term in Eq. (p6|) rapidly goes to zero. One thus obtains 



\GjW 



i(o| C »r 



(27) 



where Ij is the inverse participation ratio (IPR) for the 
Fock state Cj\FS). In other words, Pj — is the 

number of many-body eigenstates that overlap with the 
quasiparticle state c'j\FS). If a quasiparticle is localized 
in Fock space, Pj w 1, while in the extended case Pj 3> 1. 
Thus, for extended states Ij — ► in the thermodynamic 
limit. In the later case one would expect the transition 
from localized and extended states to be marked by an 
abrupt jump in Ij as some parameter (the excitation en- 
ergy, Sj, interaction strength, A, or the dimensionless con- 
ductance, g) is varied. 



III. NUMERICAL SIMULATIONS 



Equations QI9|) and (|2l|) have obvious advantages over 
Eq. (|l5|). Knowing gj(t) and fj(t) for all j and at all 
times t < to allows one to determine gj(to + St) and 
fj(to+St) by iteration, for sufficiently small St. Thus, one 
can calculate the Green's function at any time starting 
from the initial conditions gj(0) = fj(0) = 1 and gj(0) = 
fj(0) = (the latter follow from the integro-diffcrcntial 
equations because the kernels are non-singular). Since 
no inversion or diagonalization is involved, one can treat 
rather large systems. In the present method the time 
cost increases with the number M of single-particle basis 
states. For the unconstrained summations in Eqs. (|2C 



and (22) the time it takes to go through K iterations 
goes as K 2 M A when N < M. Nevertheless, this time 



can be substantially shorter than the diagonalization of 
a full many-body Hamiltonian matrix defined over the 
same M single-particle basis states. 

The weakest point of the iterative calculation is the 
appearance of spurious numerical instabilities caused by 
the non-linearity of the coupled equations. However, we 
have found a good balance between K, N, and g (for sim- 
plicity, we fixed A = 1). When the interactions are strong 
(small g) and instabilities are more pronounced, the func- 
tions gj(t) and fj(t) decay faster and, despite St being 
small (as compared to the Heisenberg time tn = 1/A), 
K is not exceedingly large. If the interactions are weak 
(large g), the decay is slow; we thus need to keep iterat- 
ing up to t tn to achieve the stationary state implied 
by Eq. (p7|). However, in this case instabilities nearly do 
not occur and we can use larger St. The number K was 
thus approximately the same for all values of g. 

As pointed out in Ref. [l3|, one can simplify the nu- 
merical calculations without sacrificing the physics by 
restricting the number of connected states in Fock space. 
In the so-called layer model,Eil all non-diagonal matrix 
elements Vj™ except those with j + k = I + m are set to 
zero. This model provides a good approximation when 
the typical non-zero matrix element is much smaller than 
A. For the layer model the number of total operations 
after K iterations decreases to K 2 M 3 . 

In the numerical simulations we also assume that the 
single-particle energy levels are random variables uni- 
formly distributed in a sequence of intervals 



Ej e [A(i-l/2);A(i + l/2)] 



(28) 



j = 1,2,.. .A/Jii3 This has the clear advantage of pro- 
viding some degree of level repulsion, consistent with 
Gaussian orthogonal (time-reversal symmetric) ensemble 
of random matrices, but without the inconvenience of a 
non-uniform level density. The number N of states with 
i < defines the number of particles present in the Fermi 
sea. The matrix elements of the interaction have to obey 
certain symmetries, namely, 



vir = V kl 

kl jrn 



VP" 

v kl 



-v k T 



-V, J " 

v lk 



(hcrmiticity) 
(fermionic character). 



(29) 



Each independent matrix element is chosen from a Gaus- 
sian distribution with zero average and SV = A/g stan- 
dard deviation. Notice that, after setting A = 1, the 
interaction strength is controlled by the dimensionless 
conductance g > 1. The numerical algorithm solves the 
integro-differential equation by iteration (Cauchy prob- 
lem) , as in simulations of the time evolution of dynamical 
systems. 



A. Results 

Below we present the numerical results obtained for 
the layer model with 50 particles and 150 basis state and 
various values of g. 
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Figure || shows the averaged return probability | Gj (t) \ 
as a function of time t for various initial single-particle 
excited states j when g = 100. For comparison, we have 
also plotted the same quantity taken for a single real- 
ization (no averaging). Notice that in the latter case 
the return probability decays rapidly and then oscillates 
strongly. The averaging over different realizations sup- 
presses such oscillations at large times, yielding a satu- 
ration at the inverse participation ratio value Ij. The 
higher the initial state energy Sj, the lower the inverse 
participation ratio, showing an increasing coupling to the 
many-body eigenstates. 




1500 



FIG. 2. Decay of the average return probability | Gj (i ) | 2 for 
g = 100 and various values of Ej. The solid line corresponds 
to a single realization for j = 20. Time is given in units of 
Heisenberg time, tH- 

Figure || shows the local density of states (LDOS) 
Pj(E) — —^lmGj(E) for a given realization of the sys- 
tem, but different values of j. This quantity character- 
izes the spreading of the single-particles states over the 
many-body eigenstates. One may observe three distinct 
regimes as the initial state energy is increased. At low 
j the initial state couples only to one or a few many- 
body eigenstates, thus corresponding to a state localized 
in Fock space. As the energy Sj increases, the initial 
state begins to couple to a larger number of many-body 
eigenstates; the resulting LDOS shows a "fractal" struc- 
ture, with multiple peaks located near zero energy. As 
Sj increases further, the initial single-particle state over- 
laps to a much large number of many-body states, result- 
ing in the merging of peaks and a rather broad LDOS. 
It is important to remark that upon ensemble averag- 
ing, pj (E) always develops into a Breit-Wigner distribu- 
tion (Lorentzian) regardless of the value of j. Only the 
width of the distribution depends on the initial state en- 
ergy. Using the GR, one can estimate this width to be 

According to the CT model predictions the transi- 
tion between localized and extended states in Fock space 
occurs at energy E* « Ayg/ lng. This implies a dereal- 
ization threshold at E*/A w 4.7 for g = 100. However, 



Fig. |H suggests instead that quasiparticles survive for 
values of Sj up to 10. In order to determine the posi- 
tion of the derealization threshold, we have adopted the 
following criterion. A single-particle state is considered 
well defined when its coupling to many-body states is 
very weak and the corresponding participation number 
Pj should not exceed some small value (of the order of 
unity). Thus, we have estimated E* as the energy neces- 
sary for Pj to reach a fixed value (for instance, 2). The 
resulting scaling behavior with respect to g is shown in 
Fig. U and is accurately described by E* ~ A ^fg. We also 
obtain the same scaling fixing Pj — 10. The difference 
between the square root scaling and the CT prediction is 
likely to be linked to the assumption of constant coordi- 
nation number made in Ref. |H[ As in the real situation, 
the branching number at each generation varies in the 
numerical simulations; eventually, it decreases abruptly 
when the number of particle-hole states is of the order 
of the number of particles in the system. As it was sug- 
gested in Ref. |[ taking into account the decreasing of the 
branching number with increasing the generation shifts 
the threshold energy to higher levels. This agrees with 
our numerical data. 
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FIG. 3. Local density of states Pj(E) = —-TmGj(E) for 
g — 100 but different quasiparticles energies Ej. F is the 
Breit-Wigner width expected from the golden rule prediction. 

Following Ref. |l2|, we attempt to characterize the point 
of transition between the different regimes by plotting 
the average participation number Pj as a function of j 
(energy) for different values of g. The result is shown 
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in Fig. H The absence of cusps or shoulders in the 
graphs is a clear indication that no sharp transition be- 
tween localized and extended states in Fock space takes 
place. Instead, the data is consistent with a continu- 
ous crossover, as previous exact diagonalization of small 
systems indicated.Ej We have checked that the same be- 
havior holds for systems twice as large as well. 



100 



small- rr. perturbation theory diverges and F is expected 
to follow a behavior that can be deduced from the GR. 




10000 



FIG. 4. Derealization energy E* for which Pj = 2 (circles) 
or Pj — 10 (squares) as a function of g. The dotted lines show 
the fits for E* ~ ^fg. 

In the fully delocalized, ergodic regime, quasiparticles 
couple to all energetically allowed many-body states. A 
way to characterize this regime is to fit the IPR Ij by 
the GR prediction, namely, Ij « min(l,<5/r), where S 
is the mean maw-body level spacing at a given energy. 
Since S/T oc g 2 fS one expects a quadratic dependence of 
Ij with g. In Fig. ^| we present the IPR Ij as function 
of g for different values of j. It is clearly visible that 
there is a faster than quadratic increase in the IPR for 
high values of j. The quadratic behavior is restored as 
the interaction strength grows (g decreases). The failure 
of the GR prediction at large g occurs as one crosses a 
certain critical value, putting a given energy Sj into the 
non-ergodic regime, namely, Sj < E**(g). Although less 
pronounced, a similar behavior was observed for smaller 



systems in Rcf. 13 



While Fig. H yields some information about the region 
where the IPR deviates from the GR prediction, it does 
not make clear where exactly ergodicity is broken. To 
quantify that Leyronas et alUB introduced a rescaled form 
of the IPR, F = —ylnlj, where x = (e/gA)lng is the 
rescaled dimensionless-conductance and y — g(A/e) 3 ^ 2 
is the rescaling factor .EJ This quantity has the advantage 
of illustrating both extreme situations of localized and 
ergodic (delocalized) regimes in a much clearer way than 
Ij itself. For x < 1 states are localized in Fock space and 
one can expand F in powers of x for a fixed Ej ; Rcf. |lj 
argues that the transition to the ergodic regime occurs 
for values of the parameter x ~ 1. In fact, at x ^> 1, the 
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FIG. 5. Participation ratio as function of energy Ej 
ferent values of g. 
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FIG. 6. Inverse participation ratio Ij as function of the di- 
mensionless conductance g for different values of the energy 
Ej. The dotted lines indicate the Fermi golden rule prediction. 

In Fig. [?] we have plotted F as function of the rescaled 
dimensionless conductance x. Our results agree with 
those presented in Ref. [15] for small energies. However, 
as the energy is increased, three regimes (instead of two, 
as seen in Ref. [l5]) clearly appear. For large values of g 
or small x (x < 1), where the perturbative calculation 
holds, the maximum of the function F shifts towards 
x = 0, as Sj is increased. This may be interpreted as 
the appearance of higher than linear contributions to the 
perturbative expansion. Another notable behavior of the 
graphs at high Ej is the absence of matching to the GR 
prediction for an extended range of x (when x > 1). For 
instance, when j = 90, the GR prediction agrees with 
the numerical data only at x > 4. Thus, our results 
indicate existence of a range of values of g where there 
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is neither localization nor ergodicity (full derealization) . 
This range grows in size at higher energies. 




x=(£/gA)lng 

FIG. 7. Rescaled inverse participation ratio as function of 
the effective interaction strength x at different energies. The 
dashed lines are the golder rule prediction. The inset shows 
in the detail the large- a; behavior. 



IV. DISCUSSION AND CONCLUSIONS 

Using an approximate method in the time domain we 
have calculated the single-particle Green's function and 
extracted information about the quasiparticle decay rate. 
The method allows us to study larger systems than previ- 
ously treated in the literature. At low energies our results 
are comparable to those obtained from direct diagonal- 
ization in small systems. At higher energies, however, 
we are able to identify three different regimes as function 
of the effective interaction strength. Furthermore, we 
have found that the derealization threshold occurs at 
energies higher than the prediction of the CT mapping, 
Ay^g/ In g; from our data, we estimate E* ~ A^fg. 

It is interesting to notice that a sequence of three 
different regimes in the absence of sharp transitions 
also occurs for quantum systems whose classical dy- 
namics corresponds to the Komolgorov-Arnold-Moser 
(KAM) regimeO Indeed, drawing an analogy between 
the present problem and the KAM theorem of classical 
mechanics is very tempting. In that sense, the prob- 
lem of a quasiparticle decay in a finite interacting system 
would be interpreted as a type of soft, interaction-induced 
crossover to chaosll3 rather than an abrupt transition for 
a particular value of £j . 

Within this interpretation, the scaling behavior seen in. 
the perturbative expansion around the localized stateEl 
corresponds to the theory of a regular classical sys- 
tem subjected to small symmetry breaking (potentially 
chaotic) perturbation. The radius of convergence of such 
theory breaks down at finite perturbation strength, indi- 
cating the onset of mixed dynamics (destruction of the 



first phase space tori). Beyond the convergence radius, 
soft chaos is present, but there is no complete ergodicity 
yet. Indeed, our data for the IPR follow a similar be- 
havior for intermediate values of the interaction strength 
g. The KAM analogy would be completed by connect- 
ing the onset of ergodicity in classical phase space (at 
sufficiently strong perturbations) to onset of the regime 
where the GR prediction for the IPR apply. 

However, it is important to remark that we have not 
attempted to construct in any form a classical counter- 
part of the quantum many-body system.E2l Also, our ap- 
proach does not give any information about the nature 
or the statistical fluctuations of the many-body excited 
states. Thus, we cannot confirm whether the quasipar- 
ticle decay problem and the KAM theorem are deeply 
related, or they just provide coincident regimes. 

Recently, a few theoretical predictions have been 
made ajbruit the dynamics that governs the decaying 
processEj'EJ So far, no experimental or numerical simula- 
tions have corroborate these predictions. We plan to in- 
vestigate this problem using the method here introduced. 
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